% This function computes the nominal and real term premia. For 
% TP_version == 1 we use the defition of term premium in Rudebusch and
% Swanson. For TP_version = 2 we use the definition in Dai and Singleton
% (2002) where the convexity term is included in the term premium

function out = getTermPremia_loadings(model,setupFilter,TP_version)

Max_maturity = max(setupFilter.maturities); %maximum maturity
ne           = size(model.eta,2);
nx           = size(model.hx,1);
ny           = size(model.gx,1);
vectorMom3   = zeros(1,ne);

%% The Nominal Term Premium
if TP_version == 1
    % Approximation of the bond prices under Q
    %The position of the first bond price in the g-function for modelFull
    ny_1bond     = find(strcmp(model.labely,'$p1_t$')); 
    nM           = exp(model.g0(ny_1bond));
    nMxp         = zeros(1,nx);
    nMyp         = zeros(1,ny);
    nMxpx        = zeros(nx,nx);
    nMxpxp       = zeros(nx,nx);
    nMxpy        = zeros(nx,ny);
    nMxpyp       = zeros(nx,ny);
    nMypx        = zeros(ny,nx);
    nMypxp       = zeros(ny,nx);
    nMypy        = zeros(ny,ny);
    nMypyp       = zeros(ny,ny);
    if setupFilter.MEXon == 1
        [Q,Qx,Qxx,Qxxx,Qss,Qssx,Qsss] = Get_Bond_Prices_Mom3_LogMex(nM,nMxp,nMyp,...
            nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
            model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
            model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,model.eta,...
            vectorMom3,int32(Max_maturity),int32(ny_1bond),int32(setupFilter.appOrder),...
            int32(1),int32(ny),int32(nx),int32(ne));
    else
        [Q,Qx,Qxx,Qxxx,Qss,Qssx,Qsss] = Get_Bond_Prices_Mom3_Log(nM,nMxp,nMyp,...
           nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
           model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
           model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
           model.eta,vectorMom3,Max_maturity,ny_1bond,setupFilter.appOrder,1);
    end        
    % Computing the nominal term premia
    % Allocating memory for the output
    TPx   = zeros(Max_maturity,nx);
    TPxx  = zeros(Max_maturity,nx*nx);
    TPss  = zeros(Max_maturity,1);
    TPxxx = zeros(Max_maturity,nx*nx*nx);
    TPssx = zeros(Max_maturity,nx);
    TPsss = zeros(Max_maturity,1);
    % Term premia = y_P - y_Q = -1/n*logP^n_t - (-1/n*logQ^n_t)
    %                         = -1/n*(logP^n_t - logQ^n_t)
    %                         = -1/n*(-n/1*y_P^(n)_t - logQ^n_t)
    for i=1:Max_maturity
        TPx(i,:)       = -1/i*(-i/1*model.byx(i,:)-Qx(i,:));
        TPxx(i,:)      = -1/i*(-i/1*model.byxx(i,:)-reshape(Qxx(i,:,:),1,nx^2));
        TPss(i,1)      = -1/i*(-i/1*model.byss(i,1)-Qss(i,1));
        TPxxx(i,:)     = -1/i*(-i/1*model.byxxx(i,:)-reshape(Qxxx(i,:,:,:),1,nx^3));
        TPssx(i,:)     = -1/i*(-i/1*model.byssx(i,:)-Qssx(i,:));
        TPsss(i,1)     = -1/i*(-i/1*model.bysss(i,1)-Qsss(i,1));
    end
else
    % Computing the conditional expected short rates
    y_select = find(strcmp(model.labely,'$r_t$')); % The short rate is stored at position at element 2
    % The arrays contain expectations from 1,2,..., Max_maturity into the future
    [rx,rxx,rss,rxxx,rssx,rsss] = CondMoments_Mom3_log(model.gx,model.gxx,model.gss,...
        model.gxxx,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,y_select,Max_maturity-1,setupFilter.appOrder);

    % Computing the nominal term premia
    % Allocating memory for the output
    TPx   = zeros(Max_maturity,nx);
    TPxx  = zeros(Max_maturity,nx*nx);
    TPss  = zeros(Max_maturity,1);
    TPxxx = zeros(Max_maturity,nx*nx*nx);
    TPssx = zeros(Max_maturity,nx);
    TPsss = zeros(Max_maturity,1);
    % Term premia = R_long - E_t*average(R_t+i-1,i=1:n))
    % For the nominal term premia
    cum_sumx           = zeros(Max_maturity,nx);
    cum_sumxx          = zeros(Max_maturity,nx*nx);
    cum_sumss          = zeros(Max_maturity,1);
    cum_sumxxx         = zeros(Max_maturity,nx*nx*nx);
    cum_sumssx         = zeros(Max_maturity,nx);
    cum_sumsss         = zeros(Max_maturity,1);
    % The first term is special:
    cum_sumx(1,:)      = model.gx(y_select,:);
    cum_sumxx(1,:)     = reshape(model.gxx(y_select,:,:),1,nx^2);
    cum_sumss(1,1)     = model.gss(y_select,1);
    cum_sumxxx(1,:)    = reshape(model.gxxx(y_select,:,:,:),1,nx^3);
    cum_sumssx(1,:)    = model.gssx(y_select,:);
    cum_sumsss(1,1)    = model.gsss(y_select,1);
    for i=2:Max_maturity
        cum_sumx(i,:)       = cum_sumx(i-1,:)    + rx(i-1,:);
        cum_sumxx(i,:)      = cum_sumxx(i-1,:)   + reshape(rxx(i-1,:,:),1,nx^2);
        cum_sumss(i,1)      = cum_sumss(i-1,1)   + rss(i-1,1);
        cum_sumxxx(i,:)     = cum_sumxxx(i-1,:)  + reshape(rxxx(i-1,:,:,:),1,nx^3);    
        cum_sumssx(i,:)     = cum_sumssx(i-1,:)  + rssx(i-1,:);
        cum_sumsss(i,1)     = cum_sumsss(i-1,1)  + rsss(i-1,1);    
    end
    for i=1:Max_maturity
        TPx(i,:)       = model.byx(i,:)       - 1/i*cum_sumx(i,:);
        TPxx(i,:)      = model.byxx(i,:)      - 1/i*cum_sumxx(i,:);
        TPss(i,1)      = model.byss(i,1)      - 1/i*cum_sumss(i,1);
        TPxxx(i,:)     = model.byxxx(i,:)     - 1/i*cum_sumxxx(i,:);
        TPssx(i,:)     = model.byssx(i,:)     - 1/i*cum_sumssx(i,:);
        TPsss(i,1)     = model.bysss(i,1)     - 1/i*cum_sumsss(i,1);    
    end
end

%% The Real Term Premium
% Solve for real bond prices
[nM,nMxp,nMyp,nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp] = numDerivM_Real(model.params);
vectorMom3   = zeros(1,ne);
ny_1bond     = find(strcmp(model.labely,'$p1Real_t$'));
solveRiskTermsSmallModel = 1;
if setupFilter.MEXon == 1
    [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3_LogMex(nM,nMxp,nMyp,...
    nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
    model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
    model.eta,vectorMom3,int32(Max_maturity),int32(ny_1bond),int32(model.appOrder),int32(solveRiskTermsSmallModel),int32(ny),int32(nx),int32(ne));
else
    [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3_Log(nM,nMxp,nMyp,...
        nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
        model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
        model.eta,vectorMom3,Max_maturity,ny_1bond,model.appOrder,solveRiskTermsSmallModel);
end
% Getting Real Yields
model.byReal0    = -bsxfun(@ldivide,(1:Max_maturity)',log(P(:,:))); 
model.byRealx    = -bsxfun(@ldivide,(1:Max_maturity)',Px(:,:));
model.byRealxx   = -bsxfun(@ldivide,(1:Max_maturity)',Pxx(:,:,:));
model.ByRealxx   = -0.5*reshape(bsxfun(@ldivide,(1:Max_maturity)',Pxx(:,:,:)),Max_maturity,nx^2);
model.byRealss   = -bsxfun(@ldivide,(1:Max_maturity)',Pss(:));
model.byRealxxx  = -bsxfun(@ldivide,(1:Max_maturity)',Pxxx(:,:,:,:));
model.ByRealxxx  = -1/6*reshape(bsxfun(@ldivide,(1:Max_maturity)',Pxxx(:,:,:,:)),Max_maturity,nx^3);
model.byRealssx  = -bsxfun(@ldivide,(1:Max_maturity)',Pssx(:,:));
model.byRealsss  = -bsxfun(@ldivide,(1:Max_maturity)',Psss(:));

if TP_version == 1
    % Approximation of the bond prices under Q
    %The position of the first bond price in the g-function for modelFull
    ny_1bond     = find(strcmp(model.labely,'$p1Real_t$')); 
    nM           = exp(model.g0(ny_1bond));
    nMxp         = zeros(1,nx);
    nMyp         = zeros(1,ny);
    nMxpx        = zeros(nx,nx);
    nMxpxp       = zeros(nx,nx);
    nMxpy        = zeros(nx,ny);
    nMxpyp       = zeros(nx,ny);
    nMypx        = zeros(ny,nx);
    nMypxp       = zeros(ny,nx);
    nMypy        = zeros(ny,ny);
    nMypyp       = zeros(ny,ny);
    if setupFilter.MEXon == 1
        [Q,Qx,Qxx,Qxxx,Qss,Qssx,Qsss] = Get_Bond_Prices_Mom3_LogMex(nM,nMxp,nMyp,...
            nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
            model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
            model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,model.eta,...
            vectorMom3,int32(Max_maturity),int32(ny_1bond),int32(setupFilter.appOrder),...
            int32(1),int32(ny),int32(nx),int32(ne));
    else
        [Q,Qx,Qxx,Qxxx,Qss,Qssx,Qsss] = Get_Bond_Prices_Mom3_Log(nM,nMxp,nMyp,...
           nMxpx,nMxpxp,nMxpy,nMxpyp,nMypx,nMypxp,nMypy,nMypyp,...
           model.gx,model.gxx,model.gxxx,model.gss,model.gssx,model.gsss,...
           model.hx,model.hxx,model.hxxx,model.hss,model.hssx,model.hsss,...
           model.eta,vectorMom3,Max_maturity,ny_1bond,setupFilter.appOrder,1);
    end        
    % Computing the nominal term premia
    % Allocating memory for the output
    RealTPx   = zeros(Max_maturity,nx);
    RealTPxx  = zeros(Max_maturity,nx*nx);
    RealTPss  = zeros(Max_maturity,1);
    RealTPxxx = zeros(Max_maturity,nx*nx*nx);
    RealTPssx = zeros(Max_maturity,nx);
    RealTPsss = zeros(Max_maturity,1);
    % Term premia = y_P - y_Q
    for i=1:Max_maturity
        RealTPx(i,:)       = -1/i*(-i/1*model.byRealx(i,:)-Qx(i,:));
        RealTPxx(i,:)      = -1/i*(-i/1*model.byRealxx(i,:)-reshape(Qxx(i,:,:),1,nx^2));
        RealTPss(i,1)      = -1/i*(-i/1*model.byRealss(i,1)-Qss(i,1));
        RealTPxxx(i,:)     = -1/i*(-i/1*model.byRealxxx(i,:)-reshape(Qxxx(i,:,:,:),1,nx^3));
        RealTPssx(i,:)     = -1/i*(-i/1*model.byRealssx(i,:)-Qssx(i,:));
        RealTPsss(i,1)     = -1/i*(-i/1*model.byRealsss(i,1)-Qsss(i,1));
    end
else
    % Computing the conditional expected short rates
    % We use the real bond price and multiply it by -1 to get the real rate.
    y_select = find(strcmp(model.labely,'$p1Real_t$')); 
    % The arrays contain expectations from 1,2,..., Max_maturity into the future
    [rx,rxx,rss,rxxx,rssx,rsss] = CondMoments_Mom3_log(-model.byRealx,-model.byRealxx,-model.byRealss,...
        -model.byRealxxx,-model.byRealssx,-model.byRealsss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,y_select,Max_maturity-1,setupFilter.appOrder);

    % Computing the nominal term premia
    % Allocating memory for the output
    RealTPx   = zeros(Max_maturity,nx);
    RealTPxx  = zeros(Max_maturity,nx*nx);
    RealTPss  = zeros(Max_maturity,1);
    RealTPxxx = zeros(Max_maturity,nx*nx*nx);
    RealTPssx = zeros(Max_maturity,nx);
    RealTPsss = zeros(Max_maturity,1);
    % Term premia = R_long - E_t*average(R_t+i-1,i=1:n))
    % For the nominal term premia
    cum_sumx           = zeros(Max_maturity,nx);
    cum_sumxx          = zeros(Max_maturity,nx*nx);
    cum_sumss          = zeros(Max_maturity,1);
    cum_sumxxx         = zeros(Max_maturity,nx*nx*nx);
    cum_sumssx         = zeros(Max_maturity,nx);
    cum_sumsss         = zeros(Max_maturity,1);
    % The first term is special:
    cum_sumx(1,:)      = model.byRealx(y_select,:);
    cum_sumxx(1,:)     = reshape(model.byRealxx(y_select,:),1,nx^2);
    cum_sumss(1,1)     = model.byRealss(y_select,1);
    cum_sumxxx(1,:)    = reshape(model.byRealxxx(y_select,:,:,:),1,nx^3);
    cum_sumssx(1,:)    = model.byRealssx(y_select,:);
    cum_sumsss(1,1)    = model.byRealsss(y_select,1);
    for i=2:Max_maturity
        cum_sumx(i,:)       = cum_sumx(i-1,:)    + rx(i-1,:);
        cum_sumxx(i,:)      = cum_sumxx(i-1,:)   + reshape(rxx(i-1,:,:),1,nx^2);
        cum_sumss(i,1)      = cum_sumss(i-1,1)   + rss(i-1,1);
        cum_sumxxx(i,:)     = cum_sumxxx(i-1,:)  + reshape(rxxx(i-1,:,:,:),1,nx^3);    
        cum_sumssx(i,:)     = cum_sumssx(i-1,:)  + rssx(i-1,:);
        cum_sumsss(i,1)     = cum_sumsss(i-1,1)     + rsss(i-1,1);    
    end
    for i=1:Max_maturity
        RealTPx(i,:)       = model.byRealx(i,:)       - 1/i*cum_sumx(i,:);
        RealTPxx(i,:)      = model.byRealxx(i,:)      - 1/i*cum_sumxx(i,:);
        RealTPss(i,1)      = model.byRealss(i,1)      - 1/i*cum_sumss(i,1);
        RealTPxxx(i,:)     = model.byRealxxx(i,:)     - 1/i*cum_sumxxx(i,:);
        RealTPssx(i,:)     = model.byRealssx(i,:)     - 1/i*cum_sumssx(i,:);
        RealTPsss(i,1)     = model.byRealsss(i,1)     - 1/i*cum_sumsss(i,1);    
    end
end

%% Computing the Inflation Risk Premium
InflTPss  = TPss-RealTPss;
InflTPssx = TPssx-RealTPssx;

%% Saving output
out              = model;
out.TPss         = TPss;
out.TPssx        = TPssx;
out.RealTPss     = RealTPss;
out.RealTPssx    = RealTPssx;
out.InflTPss     = InflTPss;
out.InflTPssx    = InflTPssx;
out.byReal0      = model.byReal0;
out.byRealx      = model.byRealx;
out.byRealxx     = model.byRealxx;
out.ByRealxx     = model.ByRealxx;
out.byRealxxx    = model.byRealxxx;
out.ByRealxxx    = model.ByRealxxx;
out.byRealss     = model.byRealss;
out.byRealssx    = model.byRealssx;

end